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ABSTRACT 

We address the cosmological evolution of violent gravitational instability in high-redshift, 
massive, star-forming galactic discs. To this aim, we integrate in time the equations of mass 
and energy conservation under self-regulated instability of a two-component disc of gas and 
stars. The disc is assumed to be continuously fed by cold gas at the average cosmological 
rate. The gas forms stars and is partly driven away by stellar feedback. The gas and stars flow 
inward through the disc to a central bulge due to torques that drive angular momentum out- 
wards. The gravitational energy released by the mass inflow down the gravitational potential 
gradient drives the disc turbulence that maintains the disc unstable with a Toomre instability 
parameter Q ~ 1, compensating for the dissipative losses of the gas turbulence and raising 
the stellar velocity dispersion. We follow the velocity dispersion of stars and gas as they heat 
and cool respectively and search for disc "stabilization", to be marked by a low gas velocity 
dispersion comparable to the speed of sound ~ lOkms -1 . We vary the model parameters 
that characterize the accreted gas fraction, turbulence dissipation rate, star-formation rate, and 
stellar feedback. We find that as long as the gas input roughly follows the average cosmolog- 
ical rate, the disc instability is a robust phenomenon at high redshift till z ~ 1, driven by the 
high surface density and high gas fraction due to the intense cosmological accretion. For a 
broad range of model parameter values, the discs tend to "stabilize" at z ~ — 0.5 as they 
become dominated by hot stars. When the model parameters are pushed to extreme values, 
the discs may stabilize as early as z ~ 2, with the gas loss by strong outflows serving as the 
dominant stabilizing factor. 

Key words: galaxies: evolution - galaxies: formation - galaxies: haloes - galaxies: spiral - 
galaxies: star formation - methods: analytical 



1 INTRODUCTION 

According to our current understanding, high-redshift massive 
galaxies form in virialized dark-matter haloes that reside at the 
nodes of the cosmic web. When baryons are accreted into a halo 
along the filaments of this web (Dekel et al. 2009), angular- 
momentum conservation implies that most of the baryons settle into 
a disc, rotating with a circular velocity that is roughly comparable 
to the virial velocity of the halo (Fall & Efstathiou 1980; Mo et 
al. 1998; Bullock et al 2001; but see Danovich et al. 2011). If the 
disc is gravitationally unstable, with a Toomre parameter Q ~ 1 
(Toomre 1964, see below), the velocity dispersion a can be esti- 
mated from the circular velocity and the mass fraction in the cold 
disc component (e.g. Dekel, Sari & Ceverino 2009). 

Pioneering observations of massive galaxies at z ~ 2 have 
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revealed striking differences between local and high-redshift star- 
forming discs of comparable masses and sizes (Elmegreen & 
Elmegreen 2005; Genzel et al. 2008). In particular, the high- 
redshift discs of ~ 10 11 Mq extending to ~ 10 kpc are thick, per- 
turbed and highly turbulent, with a ~ 30—80 km s _1 , compared to 
the thin, more uniform local discs with a ~ lOkms" 1 . The high- 
redshift discs are gravitationally unstable, showing large transient 
perturbations and bound clumps of ~ lO 9 A/0 and sizes ~ 1 kpc, 
as opposed to the more uniform mass distribution in typical local 
discs, in which the molecular clouds are smaller by a few orders of 
magnitude. While all discs may be gravitationally unstable to some 
degree, the instability in the high-redshift discs is characterized by 
higher velocity dispersion and more massive perturbations and thus 
by a more rapid dynamical evolution, such as inward migration on 
an orbital timescale, which we term "violent gravitational disc in- 
stability". This is as opposed to the secular evolution associated 
with spiral arms and bars in low-redshift discs. 

The violent instability of high redshift discs can be explained 
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in the framework of a steady intense gas supply from the cosmic 
web and a rapid mass migration due to gravitational instability in 
the disc. The gravitational fragmentation of gas-rich, thick, turbu- 
lent discs into clumps and the subsequent migration into a central 
bulge have been proposed (van den Bergh et al. 1996; Elmegreen & 
Elmegreen 2005; Genzel et al. 2008; Bournaud et al. 2008) and suc- 
cessfully simulated for idealized discs in isolation (Noguchi 1999; 
Immeli et al. 2004; Bournaud et al. 2007). Violent gravitational 
disc instability in the full cosmological context has been simulated 
by Agertz et al. (2008) and Ceverino, Dekel & Bournaud (2010), 
demonstrating self-regulated instability in steady state for several 
Gyr. Dekel, Sari & Ceverino (2009, hereafter DSC09) applied a 
Toomre instability analysis (Toomre 1964) to high-redshift discs 
under the assumption that they are made of one cold component 
and fed at the average cosmological accretion rate. This analysis led 
to the tentative conclusion that an unstable disc at high redshift re- 
mains unstable in a cosmological steady state. This fails to account 
for the evolution from violently unstable discs at high redshift to 
marginally unstable or stable discs at low redshift, as indicated by 
observations and numerical simulations (e.g. § 9 in Ceverino et al. 
2011; Martig et al. 2009). 

In this paper, we generalize the DSC09 analysis by allowing 
the gas to continuously form stars. We apply a two-component disc 
instability analysis (Jog & Solomon 1984; Rafikov 2001; Wang & 
Silk 1994; Romeo & Wiegert 2011). The disc is assumed to be 
fed by fresh gas at a constant fraction of the average cosmological 
accretion rate and the gravitational instability in the disc induces 
mass inflow to a central bulge. The gravitational energy released 
by this inflow is assumed to drive velocity dispersion in the two 
components of the disc. With time, the dissipationless stellar com- 
ponent acquires high velocity dispersion, whereas the gas turbu- 
lence dissipates its energy on a dynamical tiumescale. We solve 
for the velocity dispersions of the two components under conserva- 
tion of mass and energy and the assumption that the disc instability 
is self-regulated at a Q ~ 1 state. At late times, as the accretion 
rate decreases and the disc becomes dominated by the hot stellar 
component, the gas is required to have a lower velocity dispersion 
in order to maintain the instability. When the required gas veloc- 
ity dispersion becomes comparable to the thermal speed of sound, 
c s ~ lOkms -1 , the pressure cannot keep decreasing to the level 
where gravitational instability is possible, and we associate it with 
the end of the instability phase. We examine whether the domi- 
nance of the stellar component and the energy balance associated 
with the dissipation of gas turbulence and inflow down the potential 
gradient can lead to stabilization before z = 0. 

This paper is organized as follows. In § [2] we introduce the 
ingredients of the model that describes the evolution of a self- 
regulated disc instability in a cosmological context. In § 3 we re- 
visit the one-component case, separately for stars or gas, and intro- 
duce the fiducial model for a two-component disc of gas and stars. 
In § 4 we address the two-componet evolution for values of the 
model parameters in a plausible range. In § 5 we discuss our results 
and draw conclusions. 

Throughout the paper, we assume a flat ACDM cosmol- 
ogy specified by the cosmological parameters (fi m , erg, n, h) = 
(0.27, 0.81, 0.96, 0.70), motivated by the WMAP7+BAO+H0 re- 
sults (Komatsu et al. 201 1). The Hubble time in Gyr at redshift z is 
denoted by tn(z). 



2 THE MODEL 

We integrate in time the equations of mass and energy conserva- 
tion under self-regulated, two-component disc instability. Before 
embarking on the details of the model, we present its three key 
components. 

As a galactic disc of gas mass M gas evolves, it obeys a simple 
mass conservation equation of the form 

Afgas — ^source Af s j n k : (1) 

where Af SO urce is a given gas accretion rate of cosmological ori- 
gin (i )2.1.1| >, and the sink term consists of gas conversion into stars, 
gas inflow into the disc centre and possible mass expulsions in out- 
flows. We assume that the sink term has the general form 

M S mk = M gas /r , (2) 

where r is a characteristic timescale. If A/ sourcc and r vary suffi- 
ciently slowly in time, then the equation converges on a timescale r 
(oc e~ l / T ) to a steady-state solution where A/ gas = 0, Afsourcc = 
Af s i n k, and A/g as — M SOUICC r. Analytic models based on variants 
of eq. $2% have been attempted for various purposes. Bouche et al. 
(2010) and Rrumholz & Dekel (2011) considered draining by star 
formation but ignored migration. On the other hand, DSC09 con- 
sidereed migration but ignored star formation and outflows. Dave et 
al. (201 1) considered star formation and feedback as well as the re- 
turn of outflowing gas back to the disc. Here, we consider the three 
sink terms of mass inflow towards the galaxy centre, star formation, 
and gas outflows due to stellar feedback. 

A second important component in our model is the explicit ap- 
peal to energy conservation in the context of a self-regulated disc 
instability (Sf272j. This is a new element in comparison to earlier 
studies (DSC09; Bouche et al . 2010) and it is in line with the ap- 
proach of Krumholz and Burkert (2010). At high redshift, galactic 
discs are expected to have a high surface density, reflecting the high 
mean cosmological density. The high gas accretion rate results in a 
high gas fraction, especially at t < t, when the sink term cannot yet 
catch up with the source term. Together, they lead to gravitational 
instability. In the perturbed discs, gravitational torques drive angu- 
lar momentum out and thus cause mass inflow towards the cen- 
tre, Minflow, (Gammie et al. 2002; DSC09; Krumholz and Burk- 
ert 2011; Bournaud et al. 2011). This inflow down the potential 
well between the disc outskirts and its centre, which is on the order 
of the circular velocity squared, Vj rc , provides an energy gain of 
A/inflow K^ rc , which can be used to stir up turbulence in the gas, 
characterized by a velocity dispersion <7 gas . Assuming that this en- 
ergy gain roughly compensates for the dissipative losses of the tur- 
bulence, one can write 

»> -r/2 AigaaUgas 

JWinflow Vcirc ^ . (i) 

£dis 

Here, t c u s is the timescale over which the internal turbulent energy 
of the gas is dissipated. It is expected to be comparable to one or a 
few disc dynamical times, defined as td yn = = -Rdisc/Kirc, 
where J?disc and Q are the effective disc radius and angular veloc- 
ity. The dissipation rate thus determines the disc inflow rate that is 
needed for self-regulating the disc instability at Q ~ 1. 

The third key and new element of the model is the enforce- 
ment of self-regulated instability in a two-component disc made of 
gas and stars (5f23j. A disc is unstable once the self-gravitational 
attraction acting on a patch (associated with the surface density of 
cold material E) overcomes the reppeling forces due to pressure 
(asociated with a) and rotation (associated with Q.) (Toomre 1964). 
This is expressed in terms of the Toomre parameter Q ~ crfi/E 
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being smaller than unity. If the inflow in the disc is driven by the 
instability, and in turn it provides the power necessary for stirring 
up a, the disc tends to self-regulate itself at Q ~ 1 as follows. 
When Q drops well below unity, because of a high E and/or a low 
a, the disc becomes highly prerturbed, which increases the inflow 
rate. This depletes E and stirs up a, pushing Q to above unity. As 
a result, the fragmentation stops, the inflow slows down, so E piles 
up by new accretion and a cools down, driving Q to below unity, 
and so on. The new element incorporated in our current analysis 
is treating the instability in a disc made of two components, one 
dissipative and the other dissipationless, including the continuous 
convertion of gas to stars. The generalized instability analysis de- 
scribed in ij2.3l refers to a two-component analog of the Toomre pa- 
rameter, Q2c, which is a function of the same Q, the surface densi- 
ties of the two components, and their different velocity dispersions, 
""stars > o"gas- Since <T st ars cannot decrease and may actually in- 
crease in time, the maintenance of Q 2c ~ 1 has to rely on a sas 
being low. As the gas mass fraction in the disc gradually declines, 
(jgas has to get lower and lower. However, when it becomes compa- 
rable to the thermal speed of sound, c s ~ lOkms -1 , the pressure 
cannot keep decreasing as necessary, and Q is forced to be above 
unity. We refer to this as the end of the violent instability phase, or 
"stabilization", and attempt to find out when it occurs for different 
choices of our model parameters. 

2.1 Mass Conservation 

Following eq. (TJ, the rate of change of gas mass in the disc is 
assumed to be given by 

Mgas.disc — M ga s,acc — Mgas, inflow — ^^SFr(1 + 7out) • (4) 

Here, the source term, M gas ,acc, is the gas supply rate into the disc. 
The sink term M gas , inflow is the gas mass inflow rate in the disc 
towards the central bulge, which we sometimes loosly refer to as 
"migration". The other sink term, Msfr, includes the star forma- 
tion rate (SFR), and the outflow rate due to stellar feedback, which 
is assumed to be proportional to the SFR with 7 ou t a parameter of 
order unity. Similarly, the rate of change of disc stellar mass is 



Af s , 



+ M sr 



(5) 



where M s tar,acc represents the stellar fraction in the accretion rate 
into the disc, A/ Star , inflow is the stellar mass inflow rate within the 
disc, A/sfr is the same as in eq. © but with an opposite sign, and 
the stars are assumed not to be affected by outflows. The rate of 
change of the total disc mass is the sum of 10 and {3J. The source 
terms of accretion and the SFR are estimated as described below. 
Then, by appealing to energy conservation, gas-turbulence dissi- 
pation, and self-regulated instability, we will determine the inflow 
rates within the disc, and will be in a position to integrate eq. l[4j 
and eq. (0. 

For completeness, although we do not explicitly integrate it 
with time, we note that the growth of the central bulge is described 
by 



Afblg = Mblg,mer + Mbar, inflow , 



(6) 



where Af b i g ,mcr = Af ba r — M gas ,acc — M s tar,acc is the fraction of 
the total cosmological input of baryons that come in as big mergers 
that build the bulge without joining the disc, and Afbar, inflow = 
Afgas, inflow + A/ S tar, inflow is the total baryon inflow rate within the 
disc ( i]2,l,lt . We define the bulge-to-disk ratio as 

M blg 



B/D 



-^^star,disc H - ^Wgas,disc 



(7) 



Following DSC09, we appeal to the convenient parameter 
<5disc the fraction of mass in the disc component within the charac- 
teristic radius of the disc, J?disc, 

Afstar,disc + Mgas.d ISC 



(8) 



Here the total mass M to t within fidisc includes the contributions 
of gas and stars in the disc, the stellar bulge, and the dark matter 
within i?disc- The maximum possible value of 5 is /?, the fraction 
of baryons including disc and bulge within the disc radius, 

Afbar 



<P-- 



M tl 



(9) 



The ratio of disc to total baryonic mass is then Mdisc/Mbar = 
/3~ 5disc, so that a bulge-less disc corresponds to <5disc = /3, and 
a disc-less bulge is <5di sc = 0. The bulge-to-disc ratio is B/D = 
{ft — <5disc)/5disc- F° r gas and stars separately, we have 

r Afg a s,disc n c A/star, disc 



and 



(10) 



Mtot " Mtot 

with <5disc = <5gas +<5 s tar . In order to estimate the value of /3, follow- 
ing DSC09, we crudely treat the halo as an isothermal sphere with 
a mass profile M(r) oc r, so we can write M to t — AM,; r + Mb ar , 
where A is the halo spin parameter, Af v i r is the halo virial mass, 
and Afbar is the baryonic mass. This leads to /3 ~ /bar/(/bar + A), 
where /b ar is the baryonic mass fraction within the virial radius 
(which could be lower than the universal fraction due to mass 
loss in outflows). Throughout the paper we assume A = 0.05, 
/bar = 0.075, and thus /? = 0.6. 

2.1.1 Cosmological Gas Supply 

Both analytic estimates and cosmological simulations predict that 
the cosmological baryonic input funnels into high-redshift galax- 
ies through cold streams that follow the filaments of the cosmic 
web and include merging galaxies and a smoother component (e.g., 
Birnboim & Dekel 2003; Keres et al. 2005; Dekel & Birnboim 
2006; Ocvirk et al. 2008; Dekel et al. 2009). The average baryon 
input rate is well approximated by the universal baryon fraction 
times the total cosmological input rate into galactic haloes (Dekel 
et al. 2009). We estimate the corresponding timescale for accretion 
into a halo of mass M v i r at redshift z by 

Afvir 



M v 



2.1(1 + Z )3 2 - 4 M 



0.14 o 

12 Gyr , 



(11) 



where (1 + 2)3 = (1 + z)/3 and M 12 = Af vi r/10 12 M©. This 
approximation for the average specific accretion rate has been de- 
rived by fine tuning an analytic prediction based on the EPS ap- 
proximation (Neistein et al. 2006), and it has been shown to fit well 
the halo growth rate measured in the Millennium cosmological N- 
body simulation (Neistein et al. 2008; Genel et al. 2008; Fakhouri 
et al. 2010), except that the numerical coefficient and the powers in 
eq. il lb were slightly adjusted for the WMAP7 cosmological pa- 
rameters used here. Eq. dl lb is accurate to better than 5% in the 
redshift range 0.2 < 2 < 5, and is an underestimate of ~ 10% at 
2 = and 2 = 10. 

The analytic EPS prediction in the high-redshit Einstein- 
deSitter cosmological regime is actually M oc (1 + z) 5 ^ 2 . This 
is a good approximation at 2 > 1, where the redshit is related to 
the Hubble time in Gyr as (1 + 2) ~ 6.6t~ 2 ^ 3 . For the purpose of 
a toy model that will turn out useful for back-of-the-envelope esti- 
mates, we tentatively ignore the weak mass dependence of eq. (TTJ, 
and keep the same specific accretion rate for M12 = 0.56 at 
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(1 + z) = 3. By integrating eq. WW . such a halo ends up as a 
Milky-Way halo with M12 = 2 at z — 0, which is comparable to 
the final mass of the halo used in our fiducial model below. The 
toy-model version of eq. i ll It becomes 

M ~A(l + zf /2 , A-O.CKSGyr- 1 . (12) 

This implies a simple expression for the halo and galaxy mass 
growth, 



Mae" 



(13) 



with a ~ 25.4 A ~ 0.72. A similar exponential growth has 
been found using cosmological N-body simulations (Wechsler et 
al. 2002). 

In our model, we assume that a fraction 7 acc of the average 
baryonic mass input rate actually joins the disc as gas. This is a mul- 
tiple of two factors, (1) the fraction of the baryonic input in smooth 
accretion including small mergers that join the rotating disc, as op- 
posed to more massive mergers that build the bulge, and (2) the 
fraction of gas in this smoother component, as opposed to stars that 
come in with small merging galaxies and also join the disc. We thus 
write for the specific gas accrertion rate into the disc 



Mba 



M h 



Tacc t 



1 

acc 



(14) 



We test the effects of varying 7 acc about a fiducial value of 7 acc = 
0.7. This is based on the estimate of D09 from hydro cosmological 
simulations that the average fraction of incoming mass in clumps 
that lead to mergers with a mass ratio larger than 1 : 10 is about 30%. 

We have made so far three simplifying assumptions that are 
worth mentioning. First, we limit the current analysis to the simple 
case of accretion at the average cosmological rate. In reality, the ac- 
cretion rate is varying, both among different galaxies and along the 
history of each galaxy. The effects of these variations in the accre- 
tion rate will be studied in a follow-up paper. Second, the accretion 
rate onto galaxies may be sensitive to uncertain feedback processes, 
perhaps more important at lower redshifts (e.g., van de Voort et al. 
2010). A detailed modeling of this is beyond the scope of this paper. 
Third, a fraction of the accreting mass into the disc is expected to 
be already in stars, formed earlier in small merging galaxies (e.g., 
D09; Ceverino, Dekel & Bournaud 2010). These stars could partly 
contribute to a "cold" stellar component that would participate in 
the gravitational disc instability and partly to a "hot" stellar com- 
ponent, equivalent to the bulge in terms of its effect on the disc 
instability. In our current application we do not explicitly account 
for the stars that accrete onto the cold disc, and absorb the corre- 
sponding uncertainty in the value of 7 aC c, assumed to represent the 
fraction of mass that accretes onto the disc, all in gas. 



2.1.2 Star Formation & Outflows 

As summarized in Krumholz, Dekel & McKee (2011) and refer- 
ences therein, stemming from the Kennicutt-Schnidt empirical SFR 
law (e.g., Kennicutt 1998), the star formation rate can be best mod- 
eled by a universal volumetric local star formation law. When ex- 
pressed in terms of surface densities it has the form 



tsfr 



= e s fr- 



ts 



(15) 



where Esfr is the SFR surface density, E gas is the mass surface 
density of molecular gas, and ta = [37r/(32Gp)] 1 ^ 2 is the local 



free-fall time in the star forming region, derived from the local vol- 
umetric mass density p. The dimensionless parameter e s fr is the 
SFR efficiency, i.e., the fraction of gas that is transformed into stars 
per free-fall time, and has been argued based on theory and obser- 
vations to be constant in all star forming regions, at the level of 
1 to a few percent. In high redshift discs that are self-regulated at 
a Toomre instability with Q ~ 1, the free-fall time in the giant 
clumps where stars form are comparable to the global disc crossing 
time tdyn, the inverse of the disc angular velocity £1 = Vcirc/^disc, 
which is comparable to its vertical crossing time (Krumholz, Dekel 
& McKee 2011). 

In a variant of eq. l U5h Krumholz et al. (2009, hereafter K09) 
have argued that the physics of star formation within a molecular 
cloud, and in particular the variation in tg, can be encapsulated in 
the following star formation law 



Esp 



2.6 Gyr 



v -0.33 

^gas^ 
y^O.33 
^gas.85 



if 
if 



< 1 



(16) 



where E gas ,85 = E gas /(85M0 pc -2 ). The origin of these two 
regimes is that clumps in discs with E gas ,85 > 1 are in the regime 
where they have the characteristic Toomre mass, while in galaxies 
of lower surface density the clumps collapse and fragemnt until 
they reach the critical surface density of 85Mq pc -2 that char- 
acterize low-redshift molecular clouds. The galactic discs at high 
redshift are practically always in the regime where E gas ,85 > 1 
while at low redshift they enter the other regime. In relatively mas- 
sive, metal rich discs, the gas in the disc can be assumed to be all 
molecular, while in less massive discs, especially at high redshift, 
the molecular gas fraction becomes smaller than unity and should 
be considered when deriving the molecular surface density from 
the total gas surface density (Krumholz & Dekel 201 1). 

Outflows due to stellar feedback are incorporated in eq. 10 
through the term 



Mout = 7out Msfr 



(17) 



with 7 ou t assumed to be of order unity and constant with time (see 
Table 2). This parameterization is based on theoretical and obser- 
vational results (Bouche et al. 2006; 2007; 2009; Martin & Bouche 
2009). While this takes into account the important role of feedback 
in removing gas from the disc, we assume that the contribution of 
feedback to stirring up gas turbulence on galactic scales is minor 
(Joung et al. 2009; Ostriker & Shetty 201 1; but see also Hopkins et 
al. 2011). 



2.2 Energy Conservation 

The gravitational instability in the disc is associated with torques 
that drive an outward angular-momentum flux and a correspond- 
ing mass inflow towards the disc centre (Gammie 2001; DSC09; 
Krumholz & Burkert 2010; Bournaud et al. 2011). This inflow 
involves both the gas and the stars, in the form of clump mi- 
gration as well as inflow of the smoother material in the disc 
outside the clumps. We assume that the gravitational energy re- 
leased by this mass inflow down the potential well, at a rate 
-Egas, inflow + -Estar, inflow, is deposited in the velocity dispersions 
of the two disc components. This is a "gravitational heating" ef- 
fect (Birnboim & Dekel 2008; Khochfar & Ostriker 2008). While 
the stellar component does not dissipate the acquired energy, the 
internal energy of the gas component dissipates on the turbulence 
dissipation timescale at a rate E^is- This is summarized in an en- 
ergy conservation equation for the disc internal energy, 
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-^gas. inflow ~t~ -^star, inflow 



-Edis 



(18) 



In order to solve our system of equations, we need to use an ad- 
ditional physical constraint. Here we make the strong assumption, 
that the inflow velocites of the gas and stars are the same, and that 
the energy gain is split between the two components in proportion 
to their masses. An example of stellar inflow is the clumps migra- 
tion, but the inter-clump stars are also flowing in. Such a behavior 
is indicated in zoom-in cosmological simulations (Cacciato et al., 
in preparation), but for now it should be considered a working as- 
sumption, to be fine-tuned later. This assumption allows us to split 
the energy equation into two, 



a s, inflow 



- -Edis — -ESFR 
+ £SFR , 



(19) 
(20) 



where -E gaSi j nt and J3 s tar,int represent the rate of change of internal 
energy for gas and stars, respectively, and -Esfr refers to the energy 
transfer between the two components as gas turns into stars. 

The gravitational potential difference between the outer disc 
edge and the disc centre is of order V cilc . We then rewrite the energy 
equations as 



3 d(M gas ag as ) 

2 dt 

3 d(M star q B 2 tar ) 
2 dt 



Voire Afgas, inflow - g -^g as<7 gas (^dis — T SFr)(21) 



3 Mgas <X gas 

circ Jw star. inflow ~r 

2 TSFR 



(22) 



The internal velocity is associated with a one-dimensional gas ve- 
locity dispersion <r gas . It consists of a contribution from turbulence, 
oturb, and from thermal energy, through the speed of sound c s , 
Ogas = c 2 + o"turb- We assume that c s ~ lOkms -1 , correspond- 
ing to gas cooled by atomic cooling to 10 4 K. 

The turbulence dissipation timescale in eq. J2 1 1 is 
parametrized as proportional to the disc vertical crossing 
time, which for Q ~ 1 is comparable to the crossing time in the 
disc plane, namely 



tdis — 7dis£dyn 



(23) 



The value of 7di s is expected to be about unity if the turbulence 
lengthscale £turb is comparable to the disc scale height ftdisc or 
smaller, while it could be somewhat larger, 7di s — ^turb//idisc, if 
there is turbulence on scales larger than hdisc We adopt 7di s = 3 
as our fiducial value, and study the effect of vaying this parameter 
between 1 and 10 in § 4.1. 

Three further comments are in place regarding the energy 
equations. First, since the star formation timescale is two orders of 
magnitude larger than the dynamical timescale, the dissipation term 
dominates the right hand side of eq. J2U , Thus, as long as the inter- 
nal energy of the gas is varying slowely, eq. J2U is approximated by 
eq. <[3j, where the gas inflow rate is determined by the dissipation 
rate, with only a minor correction due to the presence of the stellar 
component. Second, when the turbulence becomes very weak with 
a small aturb, the dissipation timescale becomes much longer than 
the dynamical time, and the inflow rate slows down accordingly. In 
our simple model we do not explicitly incorporate this effect, and 
instead consider the time at which (j gas becomes as small as c s as 
the time of "stabilization", which we term z sta b- Third, in eq. ( 122b . 
the change in the stellar internal energy is assumed to consist of 
two terms. The second term represents the addition of newly born 
stars with a velocity dispersion that equals the instantaneous gas 
velocity dispersion, while the first term corresponds to the gravita- 
tional "heating" of the stars by the stellar mass inflow in the disc. It 



is worth recalling that the stellar heating rate is determined by the 
assumption that the disc stars flow in with the same velocity as the 
gas. 



2.3 Two-Component Gravitational Instability 

According to the standard one-component Toomre instability anal- 
ysis (Toomre 1964), a thin rotating disc of gas or stars becomes 
unstable to axisymmetric modes once the attraction due to self- 
gravity, represented by the surface density E, overcomes both the 
centrifugal force due to rotation and the pressure force associated 
with the radial velocity dispersion a. For gas, there is an additional 
contribution from thermal pressure, but the high-redshift discs un- 
der investigation here are in a regime where the velocity dispersion 
associated with turbulence is larger than the speed of sound that 
characterizes the thermal pressure. The instability is expressed in 
terms of the condition that the Toomre parameter Q is smaller than 
a critical value Qcrit of a value about unity (Toomre 1964), 



Q 



tvGT, 



< Qcrit 



1. 



(24) 



Here k is the epicyclic frequency, related to the angular circular 
velocity fi(r) by k 2 — r dQ, 2 jdr + 4Q, 2 . For a power-law rotation 
curve Voire oc r" , this becomes k 2 — 2(1 + v)Q, 2 . We adopt here 
a flat rotation curve, as seen in cosmological simulations, namely 
v = 0. We adopt an effective value Q = Vcirc/-Rdisc, where the 
circular velocity is approximated by the virial velocity, and the disc 
radius is a constant fraction of the virial radius, i?di sc = A R v i T , 
with A = 0.05 representing the halo spin parameter. Thus, the time 
evolution of k depends only on the cosmological evolution of the 
halo virial quantities (see § 3). For a thin disc, gaseous or stellar, 
Qcrit — 1. For a thick disc, the Toomre analysis is valid as long as 
the length scale of the perturbation is larger than the disc thickness 
and smaller than the disc radius, with Q cr it ~ 0.67 (Goldreich & 
Lynden-Bell 1965). 

When the disc is made of two components, gas and stars, each 
with different E and <r, one can address the axisymmetric instabil- 
ity via a two-dimensional Toomre-like parameter Q2 C which can 
be expresses in terms of Q s t a r and Q gas , each defined by eq. d24t > 
with the E and a of the corresponding component (Jog & Solomon 
1984; Romeo 1994; Wang & Silk 1994; Rafikov 2001). The dif- 
ference in the expresion for Q2 C due to the dissipative nature of 
the gas is small when we focus attention to the most unstable scale 
(Rafikov 2001), and we ignore it here. Romeo & Wiegert (2011) 
proposed the convenient expresion that we use here, 



where 



W Q s " t L + Qgai If Qstar > Qgas 
Qs~tar + W Qjas if Qstar < Qgas : 



W - 



2o" s ta 



(25) 



(26) 



Two-component instability is characterized by Q2 C < 1. The pa- 
rameter Q2c can be thought of as a combination of Q s t ar and Q gas , 
weighted by a function of the ratio of the velocity dispersions of 
gas and stars. 

In order to account for a disc thickness, associated with a verti- 
cal velocity dispersion Oi, z for the i's component (gas or stars), the 
quantities that enter eq. {25} are Qi — TiQi, with the approxima- 
tion ~ 0.8 + 0.7(<7;, z /o-;) valid for 0.5 < o^/ffj < 1 (Romeo 
& Wiegert 2011; Romeo 1994). Setting Tj to unity corresponds to 
a thin disc in that component. Here we adopt disc thicknesses that 
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1+z 

Figure 1. Cosmological evolution of the global quantities in our fiducial 
model galaxy. Shown from top to bottom are (a) the halo virial mass M v ; r , 



(b) the virial and disc radius, i? v i r and i?dis 



XR vil , with A = 0.05, (c) 



the circualr velocity V c i rc , assumed equal to the virial velocity V v ir. and (d) 
the angular velocity Q = V c i IC / TJdisc. which at high redshift is inversely 
proportional to the cosmological time. 



are determined by the assumption of isotropic velocity dispersion 
for each component, Oi l% — Oi. As a sanity check, note that under 
this assumption Q; = 1 corresponds to Q\ w 0.67, in agreement 
with the original result for a one-component thick disc (Goldreich 
& Lynden-Bell 1965). 

Note that the two-component system can be unstable, Q2c < 
1, even when each of the components has Qi > 1. For instance, in 
the thin-disc approximation, with the gas and stars having the same 
o i and the same Sj, Q^q = 1 corresssponds to Qi — 2 for each of 
the components. 



3 MODEL PREDICTIONS 

In the current paper we focus on a halo that grows by the average 
cosmological growth rate, eq. ( lilt , into a halo comparable to the 
Milky Way halo with M12 = 2 at 2 = 0. The virial relations 
between the halo virial mass, velocity and radius are 



V 

' v 1 



GM V 

Rvir 



3M vir 
4ttR 3 ; , 



= Ap, 



(27) 



where p is the average mass density of the universe at that time. In 
the ACDM cosmology, they can be expressed as 

M12 « O.eViooA 372 w 34.27?? A' 3 , (28) 

Rvir/1 Mpc, and A is the 



where V100 = Kdr/100 km/s, Ri 
modified expansion factor, 

-1/3 



A = a 



200 0.3 V0.77 



(29) 



with a = 1/(1 + z). The cosmological time evolution of 
the density parameter fi m and the Hubble constant, h — 
77 /100 kms" 1 Mpc -1 , is given by 



f2 m (a) = 



77(a) = 77 (fi A + fi m a 



-3x1/2 



(30) 



and we use the approximation proposed by Bryan & Norman 
(1998) for the cosmological time dependence of A: 



A = 18tt - 82 Q A 



39 ni 



(31) 



The halo mass growth is shown in the top panel of Fig. Q] It roughly 
follows the approximation eq. d 1 3b . reaching M12 — 0.5 at z — 2 
and M12 — 1 at z = 1. 

The second panel of Fig.[TJshows the corresponding 7? v i r and 
7?disc = A 7? v ir with A = 0.05. The third panel of Fig.[TJshows the 
evolution of the corresponding disc circular velocity, which we as- 
sume equals the virial velocity. This velocity is growing with time 
to a flat maximum at z ~ 1.1, followed by a gradual decline to- 
wards z = 0. This can be reproduced with the mass growing as in 
the approximation eq. d 1 3 1 >. 



Vwr oc (1 + z) 1 /' 2 exp(— az/3) . 



(32) 



This indeed reaches a maximum at z max ~ 3/ (2a) — 1 ~ 1.1. 
We will see below that this governs the evolution of a in the one- 
component case, and that it also has an important effect on the evo- 
lution of the velocity dispersions in the two-component case, partly 
through the dependence of the gravitational heating on the potential 
well expressed by Vcirc- 

The bottom panel of Fig. 1 shows the angular velocity Q, 
which enters in the instability Q parameter. Recall that Q = t^™- 
In the Einstein deSitter phase where A v i r is constant, and under the 
assumption of a constant spin parameter A, tdyn is proportional to 
the cosmological tim^l 



3.1 A one-component disc 

As a first application of our model, we investigate the evolution of 
instability in a one-component disc. Isolating a form eq. {24}, we 
can write for a disc with a flat rotation curve 



7rG£disc(3c 



— ; — ^disc^circQcrit ■ 

n/2 



(33) 



Inserting this expression for a in the mass and energy equations, 
we can straightforwardly integrate them. 

Fig. E] shows the evolution of a one-component disc for dif- 
ferent values of the dissipation parameter 7di s (eq. {23}). A choice 
of 7dis ~ 00 refers to a dissipationless, all stellar disc. The finite 



1 For the standard ACDM cosmology, in its Einstein-de Sitter regime, the 
virial radius and velocity are related to the virial mass as i? v i r ~ 308 kpc 
(1 + z)- 1 M^ 3 and V vil ~ 118 km/s (1 + z) 1 / 2 M^ 3 . Thus V200 ~ 



RlOO (1 + z ) 'f/ 2 , where the quantities are in units of 200 km/s and 100 kpc. 
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Table 1. One-Component Models 
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Figure 2. Evolution of a one-component disc and the role of dissipation. 
Red solid lines refer to non-dissipative stellar discs, and blue dashed and 
dotted lines refer to gaseous discs with 7 dis = 10 and 1, respectively. Top: 
Mass fraction in the disc within -Rdisc. ^disc — V^u/Vdic - Middle: Disc 
surface density S. Bottom: Velocity dispersion a and V C \ IC for comparison. 
Lower 7di s . corresponding to a higher dissipation rate, requires a lower 
velocity dispersion for maintaining Q = 1. 



values of 7di s in the range 10-1 correspond to a pure gaseous disc 
with slow and rapid dissipation, 7di s = 10, 1 respectively. For each 
choice of 7di s we show the results for two different intitial condi- 
tions, which practically converge to a single solution by z ~ 5, in 
accordance with the exponential convergence described in fj2] 

The top panel of Fig. [2] shows the time evolution of the 
disc mass fraction, <5di sc , which in the one-component analysis is 
V2Q crit (a/Vcirc). We see that <5di S c ls roughly constant in time, 
or is varying rather slowly, indicating that the system evolves in 
a quasi-steady-state, where the disc and the bulge grow together, 
as emphasized in DCS09. In the case of a stellar disc, the disc 
fraction is high, <5disc ~ 0.4, close to the maximum possible 
value of ft = 0.6, namely a bulge-to-disc ratio of 0.5. This re- 
flects the fact that without dissipation the inflow rate in the disc 
is rather slow, eq. l[3}. For lower values of 7di s , the disc fraction 
<5disc becomes lower, refelecting the higher depletion rate by inflow 
when the dissipation is faster. With 7di s = 10, the steady-state is 
with 5disc — 0.25, namely B/D ~ 1.4, and with 7di s = 1 it is 
<5disc — 0.15, corresponding to B/D ~ 3. 





7dis 


7acc 


TSFR 7out 


Dissipation 


~ oo 


0.7 







10 


0.7 







1 


0.7 






The middle panel of Fig. [2] shows how E d j ac is gradually de- 
creasing with time. This reflects the assumptions that the galaxy 
mass follows the halo mass growth, that is close to exponential 
following eq. l |13t , and that the galaxy size is a constant fraction 
of the virial radius. Then, according to the approximate behavior 
of eq. d 1 3t and in the Einstein-deSitter regime with a constant A, 
the three-dimensional density is proportional to the virial density, 
namely it is independent of mass and follows the cosmological den- 
sity evolution. We thus have: 

M 



E oc 



fl- 



oe <5disc (1 + z) exp(— az/3) . 



(34) 



For a slowly varying <5disc, as seen in the top panel of Fig. [2] and 
for a ~ 0.72 as in eq. d!3t , the surface density of eq. ( I34t is indeed 
decreasing with time at all redshifts z < 7.4, namely throughout 
most of the relevant redshift range. 

The lower panel of Fig. [2] shows the evolution of the velocity 
dispersion in the one-component models, with V c irc shown as a ref- 
erence. Based on eq. i33l . a is determined by S/fi, or equivalently 
by 5diac V c i Ic . Since 5disc is rather constant, a is roughly propor- 
tional to Kirc, which follows the evolution of the virial velocity 
shown in Fig.Q]and approximated in eq. ((32}. This explains why a 
is rising to a flat maximum near z ~ 1 and is gradualy declining 
afterwards. Thus, the evolution of a is driven by the cosmological 
halo growth rate and the evolution of its virial velocity. More pre- 
cisely, the evolution of a is determined by the interplay between 
the rates of change of E and Q under the constraint that Q ~ 1 
(eq. J24H. Figure [3] demonstrates this, and in particular the quali- 
tative role played by the sink terms in causing the decline of a at 
late times. The top panel shows a simplified model in which mass 
accretes onto the disk in the average cosmological rate, but the sink 
terms are all turned off, <5disc = const., and there is no dissipation. 
In this case, the evolution of E and f2 are dictated solely by the cos- 
mological evolution of the virial quantities, under the assumption 
of a cosntant spin parameter and Kirc = Kd r . We note that in this 
case, in the redshift range z = 3 — 0, that the decline rate of E and 
Q are quite similar, so a maintains a rather cosntant value. The toy 
model of eq. j!3K in the Eisntein-deSitter regime, indeed predicts a 
slow evolution a oc E/fl oc (1 + z) 1//2 e~ 0,24z . The bottom panel 
of Fig. prefers to the same model but with an artificial suppression 
term (oc (1 + z) 1//2 ). In this case, at late times, E drops in time 
faster than fl, so a has to decline in order to maintain Q ~ 1, as it 
does in Fig. [2] 

For the models shown in Fig. [2] we see that self-regulated in- 
stability is maintained till z — 0, in the sense that the model has 
fgas > c s ~ lOkms -1 at all times, even for the case in which 
the dissipation rate is as high as the disc crossing time. This re- 
sult is consistent with the analysis of DSC09, who concluded that 
an unstable disc accreting at the average cosmological rate remains 
unstable throughout its cosmic history. 

We note that DSC09 estimated the migration timescale via an 
explicit calculation of the effect of encounters between giant clump, 
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Figure 3. An illustration of the way the velocity dispersion compensates 
for the rate of change of E/fi to keep Q ~ 1. 7bp: a case where the disc 
grows by cosmological accretion but all sink terms are turned off, showing 
(r oc S/f! ~ const.. Bottom: the same but with an arbitrary sink term 
added, showing a decline of a in the regime where Y./Q declines, similar 
to Fig. [2] 



100 F 



CO 
O 
w 
CD 

B 



T 



t 1 It 1 ' '! 

gas. acc acc 



t innow (DSC09) 

7 dis =3 




This Work 
- [one component] 



1 



1 



6 8 



1+z 



Figure 4. Evolution of relevant timescales, showing the cosmological ac- 
cretion timescale (black dotted line), the timescale for mass inflow in the 
disc as derived by DSC09 from clump interactions and torques (blue dashed 
line), the timescale for mass inflow in the disc according to our fiducial 
model (red solid line), and its range corresponding to variations in 7^ in 
the range (1, 10) (yellow shaded area). 



Table 2. Two-Component Models 
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and alternatively via torques in a viscous disc, without consider- 
ing the energy budget. In Fig. [4] we show that our current model, 
based on energy considerations, without specifying the mechanism 
by which this energy is transferred to velocity dispersion and to 
inflow in the disc, predicts an inflow timescale similar to the esti- 
mate in DSC09. This agreement is non-trivial given the fundamen- 
tal difference between the two calculations. We also see that the 
two calculations predict inflow timescales that are comparable to 
the accretion timescale. 



3.2 Disc of Gas and Stars: Fiducial Model 

After recovering the persistance of instability in the one-component 
models, we proceed to the two-component case. Here, the gas frac- 
tion in the disc is gradually declining with time, which may lead 
to stabilization once the disc becomes dominated by a hot stellar 
component. In this section, we describe the results for a specific set 
of model parameters, which we refer to as the fiducial model. In Sj4] 
we will address the effects of varying the model parameters about 
the fiducial values. 

Our fiducial model, as specified in Table 2, consists of 7di s = 
3 (eq. (23)), 7 acc = 0.7 (eq. (HJl), the K09 star formation law 
(eq. l|16t), and tentatively 7 out = (eq. l|17|l). Fig. [5] shows the 
time evolution of the same quantities shown in Fig. [2] namely <5disc, 
E and a, but now, and in the following figures, refering to the 
two components of gas and stars in blue and red respectively, and 
to their sum in black. We immediately note in the bottom panel 
that, unlike the always-unstable one-component models, the fidu- 
cial two-component model does stabilize shortly prior to 2 — 0, at 
z ~ 0.15, as marked by <T gas droping below c s = lOkms -1 . The 
drop is rather steep, starting at 2 ~ 1 after a pretty constant <7 gas ~ 
30kms _1 till then. It is driven by the disc turning from being gas 
dominated at early times to star dominated at z ~ 1.5, as seen in 
the top and middle panels, with £ s tar/£gas — S stSlIS /8 g3lS ~ 5 by 
z = 0.15. With the stars "hot" at a star ~ 80 - 100 kms -1 in 
the range 2 ~ 1 — 0, the gas has to cool faster than in the one- 
component case in order to compensate for its decreasing fraction 
in the disc surface density and keep Q2c = 1. 

We note several additional interesting features in the evolu- 
tion of the fiducial model compared to the one-component models. 



First, the gas velocity dispersion is at the level of a g 



30 km s 



until it starts droping at z ~ 1. At z 
Vkrc/ogas ~ 6.6, reaching V c i r c/cr gaa 



2 this corresponds to 
10 at 2 = 0.15. The 
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Figure 5. Evolution of the fiducial two-component model (Table 2). Panels 
and curves are as in Fig. [2] Shown are the quantities for the gas (blue), the 
stars (red), and their sum (black). The disc stabilizes at z sta b = 0.15. 



evolution of o- gas in the two-component fiducial model is similar 
to the evolution in the corresponding one-component gas model, 



except for the drop at z < 1. At z 



r /CTg 



2.6. Sec- 



ond, the total disc fraction within i?disc is rather constant in time, 
evolving in a quasi steady state as in the one-component gas model, 
at a comparable level of 0.2 - 0.25, i.e., B/D~ 2 - 1.4, but here 
slowly rising instead of declining. The constancy of the bulge to 
disc ratio indicates that the massive bulge is not the main reason 
for the stabilization Third, the stellar surface density in the fiducial 
model is almost constant throughout the disc evolution, contrary to 
its gradual declime in the one-component stellar model. 

Figure [6] shows the evolution of the components of Q2c, 
eq. d25t . for the fiducial model. It shows the values of Q — TQ 
for the gas and for the stars, with T ~ 1.5 under the assump- 
tion of isotropy adopted here, and the value of W. We see that 
Q gas is always smaller than Q s t a rs- This is because cr gas is sig- 
nificantly smaller than <7 s t ars , and especialy so when the surface 
density is dominated by the stars. For a similar reason W evolves 
to well below unity. As a result, Q gas is dominant in determining 
Q2c (eq. ( 1251 ), top line), namely the instability is driven by the gas 
even after the disc has become stellar dominated. 



O* 

* 

E- 




Figure 6. Evolution of the instability parameters in the fiducial two- 
component model. Top: The Q parameters, eq. {25), for gas (blue) and for 
stars (red). Bottom: The weighting function W , eq. {26). With a lower value 
of Q, the gas is the main driver of instability at all times. 



4 VARIATIONS ABOUT THE FIDUCIAL MODEL 

We now explore the effects of varying the model parameters that 
govern the turbulence dissipation rate ( ^4. U , the fraction of the 
baryonic cosmological input that joins the disc as gas ( i)4.2t . the 
star formation rate ($43}, and the outflow rate ( ij4.4t . Table 2 dis- 
plays the values of the parameters used in the different models. 



4.1 Gas Dissipation 

The energy considerations that lie at the basis of our model im- 
ply that the mass inflow rate within the disc is determined by the 
gas turbulence dissipation rate (eq. (3) and 1211 ) and that in turn 
the inflow down the potential gradient is the source of energy that 
adjusts the velocity dispersions of gas and stars at the level that 
guarantees Q2c = 1. In Fig. [7] we test the effect of varying the 
dissipation parameter 7di s (eq. l |23t ) between the extreme values 1 
and 10, bracketing the fiducial value 7di s = 3, while keeping the 
other model parameters at their fiducial values. 

We notice first that the dissipation rate has only a minor im- 
pact on the qualitative behavior, and in particular on the time of 
stabilization, which shifts from z s t a b ~ 0.1 to z s t a b ~ 0.3 when 
7dia is varied from 10 to 1, corresponding to a shift from a slow 
to a fast inflow rate. A higher dissipation rate is responsible for a 
slightly earlier stabilization. The reason for the small effect is that, 
by construction, the ratio Eg aa /E s t ar is not very sensitive to the 
inflow rate, that is driven by the dissipation rate, and this ratio is 
in fact slightly higher for faster dissipation. The transition to star 
dominance varies from z ~ 1.3 to 1.6 when 7di s is varied from 
10 to 1. Note that the surface densities of the two components are 
correlated not only because of the shared inflow and the resultant 
depletion of the two components of the disc, but also because ac- 
cording to the SFR prescription E star is growing in proportion to 
E gas . On the other hand, the value of 7di s has a non-negligible 
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Figure 7. The effects of varying 7di s (eq. fl23» in the range 1 — 10 on the 
evolution of the two-component model. Panels and curves are as in Fig. \5\ 
Dotted curves refer to the fiducial model of Fig. [5] for comparison. The 
stabilization epoch is only weakly affected, but the gas dispersion velocity 
is significantly lower when the dissipation rate is higher. 
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Figure 8. The effects of varying 7 acc (eq. dl4» in the range 0.4 — 1 on the 
evolution of the two-component model. Panels and curves are as in Fig. [3] 
and Fig. [7] The stabilization epoch is only weakly affected, but the gas 
dispersion velocity is lower when the accretion onto the disc is lower. 



effect on the value of <7 gas prior to its drop at 2 < 1, which is 
~ 40 kms -1 and ~ 20 kms" 1 for 7di s = 10 and 1, respectively. 
The corresponding values of Kirc/o" gas at z = 2 are 4.2 and 7.4 
for 7dis = 10 — 1. While the drop of <r gas may start a bit later for 
7dis = 10, the fact that it starts from a higher value makes it reach 
10 km s" 1 at a slightly later time. 

We also see, in the top panel, that a higher dissipation rate, 
inducing a higher disc depletion rate, makes the total disc fraction 
smaller, from <5disc ~ 0.25 — 0.3 to ~ 0.15 — 0.2 when 7di s is varied 
from 10 to 1. These correspond to B/D }t 1 and 2 respectively. 

4.2 Effective Accretion onto the disc 

Here, we explore the effect varying 7 acc (eq. d!4t). the fraction of 
the average cosmological accretion rate that joins the disc as gas, 
using the extreme values 7 acc = 0.4 and 1 about the fiducial value 
of 0.7. Fig.[8]displays the results, showing that the stabilization time 
shifts from z sta b — to 0.4 when 7 acc varies from 1 to 0.4. 

Higher values of 7 acc correspond to higher surface densities 
of the two components, because of the SFR dependence on E gas . 
This dictates higher values of a for the two components in order 
to maintain Q 2c = 1. Both cr star /(7 gas and E star /E gas are rather 



insensitive to 7 acc , and so are, in particular, the time of turning star 
dominated and when a gas starts dropping at z ~ 1. The case of 
higher 7 acc reaches 10 kms" 1 later because it has to drop from a 
higher value at earlier times, 35 kms -1 for 7 acc = 1 compared to 
20 kms -1 in the case of 7 acc = 0.4. 

Note also in the top panel that a higher 7 acc corresponds to a 
higher disc fraction, varying from 5disc ~ 0.17 to 0.3 when 7 acc is 
varied from 0.4 to 1. 



4.3 Star Formation 

We next study the effect of varying the star formation efficiency per 
free-fall time, e s fr (eq. H5\). in comparison with our fiducial star- 
formation law of K09 eq. U61 . Fig|9]shows the results for the cases 
with e sfr = 0.01 and 0.05. 

In order to help us interpret these results, we compare in 
Fig. 10 the star formation laws, in terms of Esfr as a function of 
E gas , for the three models studied here. We see that prior to z ~ 3, 
the Esfr in the e s f r = 0.01 model is indeed slightly higher than in 
the K09 model, and it becomes lower after this time. On the other 
hand, the e s f r = 0.05 model has a significantly higher Esfr than 
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Figure 9. The effects of varying e s f r (eq. 1151 ) in the range 0.01 — 0.05 
on the evolution of the two-component model. Panels and curves are as in 
Fig. \5\ and Fig. [7] The stabilization epoch is rather sensitive to the SFR. 
The case of very high SFR is unique in the sense that it is always stellar 
dominated, and still, it maintains instability at high redshifts, to be stabilized 
only at late times. 
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Figure 10. The SFR surface density as a function of gas surface density 
for different SFR prescriptions: the fiducial K09 prescription (black line, 
eq. {16}), and the two cases shown in Fig. [9] with e s f r = 0.01 (red line) and 
0.05 (blue line). The points refer from right to left to z = 6, 5, 4, 3, 2, 1, 
with the e s f r = 0.01 extending to z = 0. The SFR in the e s f r = 0.01 case 
is higher than the fiducial model at early times, and lower at low redshifts. 



<5disc = 0.35. The reason for the rather decline of <T ga s starting at 
z ~ 2 towards stabilization at z ~ 0.6 is subtle, as no other quan- 
tity involved shows an abrupt change near that time. One thing that 
does change though is the rate of decline of E gas , which at z ~ 2 
becomes steeper than that of Q. oc (1 + z) 3 / 2 . When this happens, 
Q2c can be maintained at unity only if <T ga s declines, and this de- 



the two other models at all times. We will see that this has an impact 
on the evolution of the instability. 

The model with e s f r = 0.01 remains unstable all the way to 
z — with an almost constant disc fraction in a steady state at 
5disc — 0.25, namely B/D ~ 1.4, not too different from the one- 
component model (§ 13. II and DSC09). The reason is that the low 
SFR at late times allows the disc to never become strongly star 
dominated — it maintains E star s/E gas < 2 over a long period of 
time, from z ~ 2 to 0. This eliminates the factor that serves as the 
main driver for stabilization in the fiducial case. On the other hand, 
at very high redshift, E star is higher than in the fiducial case. This 
is because in the high E gas regime the SFR with a constant e s f r is 
higher than with the K09 law (see Fig. 10). 

The case with very efficient SFR, e s f r = 0.05, behaves 
very differently from all other models. While it stabilizes at z ~ 
0.6, somewhat earlier than the fiducial model, the very efficient 
SFR makes this disc strongly star dominated at all times, with 
E star /E gas ~ 2.8 already at z = 7, growing to ~ 10 at later 
times. Despite the stellar dominance, the disc is unstable till after 
z ~ 1. The low gas fraction involves low dissipation and therefore 
slow inflow rate in the disc, making the disc fraction grow to above 



4.4 Outflows due to Stellar Feedback 

Figure QT] shows the effects of strong and very strong outflows via 
7out = 1 and 3 (eq. (T7J). We see that the stabilization epoch shifts 
from z sta b — 0.15 when outflows are ignored to z sta b — 0.5 and 
1 for 7 ut = 1 and 3, respectively. This means that strong outflows 
could be the strongest driver towards early stabilization among the 
physical properties explored so far. The gas removal from the disc 
due to strong outflows causes a steep decline in E gas , making it 
drop below lOOM© pc -2 already prior to z ~ 1.5 for 7 out = 
3. This translates via Q2c = 1 to a low gas velocity dispersion, 
dropping below 20 kms -1 already at z ~ 2.5, from which the way 
to stabilization at 10 km s _1 is rather short. The disc fraction in the 
extreme outflow model becomes rather low, <5di sc — 0.1, namely a 
bulge-dominated system. 

With an outflow rate comparable to the SFR, 7 ou t = 1, these 
effects are less drastic. While stabilization occurs at z s t a b ~ 0.5, 
the disc fraction is rather constant at about 0.2, and the disc at z s t a b 
has E s t a r S /E gas w 4 and <T sta rs/(T ga s ~ 10, all not that different 
from the fiducial model in which outflows were ignored. 
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Figure 11. The effects of varying 7 ou t (eq. I17H in the range — 3 on the 
evolution of the two-component model. Panels and curves are as in Fig. [5] 
and Fig. [7] Stronger outflows induce earlier stabilization, but still limited to 
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Figure 12. The effects of pushing the model parameters to extreme val- 
ues favorable of early stabilization on the evolution of the two-component 
model. Panels and curves are as in Fig. \5\ and Fig. [7j Even in this case, 
stabilization occurs only after 2 = 2. 



4.5 A Case of Extreme Stabilization 

In order to explore the robustness of the instability at high redshift, 
we present here a case in which the parameters are pushed to ex- 
treme limits within the sensible range, in an attempt to obtain early 
stabilization. This model has a very low fraction of gas accretion 
into the disc (7 gas ,acc = 0.4), a high dissipation rate (7di s = 1), 
a high outflow rate (7 out = 3), and the K09 SFR law. Figure [72l 
demonstrates that even in such an extreme case the disc remains un- 
stable until z w 2. At the time of stabilization, the surface density 
is not yet star dominated, but the low accretion rate and the strong 
outflows lead to a low total surface density, of E ~ lOOM© pc -2 
already by z ~ 2. The disc fraction is low and roughly constant at 
5 ~ 0.1, i.e., the system is bulge dominated with B/D « 5. Thus, 
the disc stabilization in this case could be interpreted as driven by 
"morphological quenching" due the massive bulge (Martig et al. 
2009) rather than by the transition to a star-dominated disc. We 
conclude that stabilization before z ~ 2 is very unlikely — it may 
be achieved only with more extreme outflows or in galaxies where 
the accretion is severely suppressed compared to the cosmological 
average. 



5 DISCUSSION AND CONCLUSIONS 

We have studied the cosmological evolution of gravitationally un- 
stable galactic discs. To this aim, we have developed an analytic 
model that evolves in time the properties of a galactic disc as a se- 
quence of quasi steady-state configurations under the assumption 
of self-regulated two-component gravitational disc instability. The 
disc is assumed to be fed by fresh gas at a constant fraction of the 
average cosmological accretion rate. The gravitational instability is 
associated with mass inflow within the disc towards a central bulge. 
In the disc, gas continuously turns into stars according to a given 
star formation law. This typically leads to a characteristic pattern 
of evolution from a gas-dominated disc at high redshifts to a star- 
dominated disc at low redshifts, which is the main driver of insta- 
bility at high redshifts and the main reason for stabilization after 
z ~ 1. 

The model is based on three key constraints, namely mass and 
energy conservation and self-regulated instability. Mass conserva- 
tion accounts for cosmolgical accretion onto the disc, mass inflow 
within the disc into a central bulge, conversion of gas into stars, 
and gas outflows. Energy conservation imposes that the gravita- 
tional energy gain due to the mass inflow down the potential gra- 
dient is driving velocity dispersion in the disc gas and stars, com- 
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pensating for the dissipation of the gas turbulence on a dynamical 
timescale. The imposed self-regulated two-component instability 
then provides a connection between the velocity dispersions of gas 
and stars, which allows us to solve the set of equations and prope- 
gate the system in time. 

While confirming that a single component disc, of gas or stars, 
is likely to maintain the instability in a steady state till z = 
(DSC09), we find here that a two-component disc, where gas con- 
tinuously turns into stars, tends to (a) be unstable at z > 1 and (b) 
stabilize at z ~ — 1. The instability at high redshift is typically 
driven by the high surface density of gas and "cold" stars, which 
is due to the high density of the universe and the high cosmolog- 
ical accretion rate. The stability at low redshift is typically driven 
by the steallar dominance, which results from the declining cosmo- 
logical accretion rate, the decline of E because of the cosmological 
expansion, the continuous convertion of gas to stars, the outflows 
of gas, and the effects of inflows in the disc on disc depletion and 
bulge growth. The low gas surface density and its fraction in the 
disc forces the gas to a low velocity dispersion in order to balance 
the high stellar velocity dispersion and keep Q2c ~ 1. 

We study the impact of varying the assumed dissipation rate, 
fraction of gas accretion onto the disc, star formation law, and out- 
flow rate. When the model parameters vary within a sensible range, 
the epoch of stabilization tend to change only in a limited range, 
typically 0.1 < z < 0.6. The inclusion of strong outflows plays 
a somewhat more significant role, and it may stabilize the disc as 
early as z ~ 1 if the mass removal rate is three times the SFR. 
An extreme model, in which the parameters are all pushed to their 
most favorable values for early stabilization, is capable of shifting 
the stabilization to an earlier epoch, but only to z ~ 2. Prior to 
z ~ 2, it is difficult to avoid violent disc instability. 

In a companion paper (Genel, Dekel & Cacciato 201 1), we ex- 
plore the alternative possibility that the disc turbulence is driven by 
the cosmological in-streaming in a non-self-regulated manner (as in 
DCS09). We find that also in this case the discs tend to evolve from 
instability to stability, and that cr/Vcirc could be in the ballpark of 
the observed values at z ~ 2, with no need to appeal to feedback 
effects. However, in this case a/Vcirc is predicted to be indepen- 
dent of the actual value of the accretion rate, thus not reproducing 
the observed decline with time. When we do impose self-regulation 
at Q ~ 1, and introduce a duty cycle for the instability, we find a 
better agreement with the observational trends. 

The current analysis is simplified in several ways that justify 
a few cautionary notes. First, the baryonic accretion is assumed to 
be at the average cosmological rate, while one expects significant 
variations in the accretion rate along the history of every galaxy and 
between different galaxies as a function of their environment. Since 
the accretion rate is the main driver of disc instability and the asso- 
ciated phenomena, these variations may have important effects on 
the evolution. The evolution of instability under realistic accretion 
histories will be addressed in a follow-up paper. 

Second, in order to solve the set of equations, we had to make 
assumptions concerning the stellar inflow rate in the disc compared 
to the gas inflow rate, and especialy how the gravitational energy 
gain by the inflow is distributed between stirring up turbulence in 
the gas and raising the velocity disopersion of the stars. Our as- 
sumptions boils down to assuming that (a) the mass inflow rate per 
unit mass, or the drift velocity, is the same for gas and stars, and (b) 
that the energy deposited per unit mass in the disc is the same for 
gas and stars. The participation of stars in the inflow within the disc 
is obvious from the fact that they dominate the giant clumps that are 
known from analytic estimates (DSC09 and references therein) and 



simulations (Ceverino et al. 201 1 and references therein) to migrate 
inward on an orbital timescale as a result of torques, inter-clump 
interactions, and dynamical friction. However, the validity of the 
assumptions adopted here for the exact behavior of the stars com- 
pared to the gas is still uncertain. We note, for example, that Forbes 
et al. (in preparation), in a study of a similar problem, are adopt- 
ing a different assumption concerning the stellar migration rate (M. 
Krumholz and J. Forbes, private communication). The inflow rate 
in the two disc components, the migration of clumps compared to 
the inflow of disc mass outside the clumps, and the evolution of ve- 
locity dispersion in the two components in the context of gravita- 
tional instability are all being studied in detail using high-resolution 
cosmological simulations (Cacciato et al. in preparation). 

Third, we focused here on massive galaxies, comparable in 
mass to the Milky Way and to the observed big clumpy discs at 
z ~ 2 (Genzel et al. 2006; 2008; Elmegreen & Elmegreen 2005). 
However, one should be aware of the possibility that the evolu- 
tion of instability may be different in less massive galaxies. On one 
hand, at a given redshift, less massive galaxies tend to have a sim- 
ilar £1 but a lower total surface density, which would require lower 
velocity dispersions for maintaining Q2c ~ 1, thus leading to ear- 
lier stabilization. On the other hand, if the less massive galaxies are 
more gas rich, as observed, perhaps due to different effects of star 
formation and feedback, then stabilization is expected later in less 
massive galaxies. There are indications for such a "downsizing" in 
instability in observations (Elmegreen & Elmegreen 2005) and in 
simulations (Bournaud et al. 2011). The mass dependence will be 
addressed in an upcoming paper. 

Despite these limitations, our analysis indicates that phases of 
violent disc instability are a solid prediction at z > 1, and that 
the typical discs are likely to stabilize prior to the present epoch, 
mostly due to the growing dominance of "hot" stars at late times. 
The predictions of our analytic model are to be compared to those 
from high-resolution hydro-cosmological simulations. 
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